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Abstract. 

Mutation is introduced into autocatalytic reaction networks. The differential equa- 
tions obtained are neither of rephcator type nor can they be transformed straight- 
way into a hnear equation. Examples of low dimensional dynamical systems — 
n = 2, 3 and 4 — are discussed and complete qualitative analysis is presented. 
Error thresholds known from simple replication-mutation kinetics with frequency 
independent replication rates occur here as well. Instead of cooperative transi- 
tions or higher order phase transitions the thresholds appear here as supercritical 
or subcritical bifurcations being analogous to first order phase transitions. 



Keywords. 
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1. Introduction 

About ten years ago Schuster and Sigmund (1983) pointed out that the same 
class of model equations applies to the description of population dynamics in differ- 
ent disciplines ranging from chemical kinetics to theoretical sociobiology and they 
coined the terms replicator dynamics for the common phenomena and repli- 
cator equations for the corresponding ODEs, respectively. These simple models 
deal exclusively with (error-free) replication in populations of constant size. We 
are concerned here with more complex dynamics which allows for mutations within 
a set of variants, all of which are considered explicitly as variables. In particular we 
study some low-dimensional examples for which analytical treatments are possi- 
ble, and present a few results for arbitrary dimensions. Low-dimensional replicator 
equations with reaction channels for error copies are for example relevant models 
for introducing mutation into 

• chemical kinetics of catalyzed replication, 

• Fisher's selection equation, 

• dynamics of coevolution, and 

• dynamics of Maynard Smith games. 

It is worth mentioning that replicator equations of dimension n are dynamically 
equivalent to Lotka-Volterra equations of dimension n — 1 (Hofbauer, 1981). 

Explicit consideration of mutation terms is highly relevant for realistic mod- 
els since they describe processes of fundamental importance in nature which are 
often neglected only because they make the handling of the model equations much 
harder, and they are commonly prohibitive for analytical solutions. An often use- 
ful compromise is shortly mentioned here: variants are not described by explicit 
variables but rather lumped together into a so-called "error-tail" as suggested in a 
stochastic analysis of error propagation by Nowak and Schuster (l989). This sim- 
plification allows to treat cases with more elaborate dynamics and autocatalytic 
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reaction networks (Andrade et a/., 1993; Nuno et a/., 1993a, b; Stadler and Nuno, 
1994). 

Replication processes, correct and incorrect, are fundamental to selection and 
evolutionary dynamics in general. In the simplest case which refers to asexual 
reproduction in constant environments replication rates are proportional to tem- 
plate concentrations and the replication rate parameters are constants. Then 
populations approach a unique asymptotically stable steady state which has been 
characterized as (molecular) quasi-species since these stationary mutant distribu- 
tions represent the genetic reservoirs of asexually replicating species as much as 
the gene pools do in ordinary biological species. The formal mathematical back- 
ground of the quasi-species concept was updated and its relation to experimental 
data obtained in test tube evolution experiments and in virology was summarized 
in a recent review (Eigen et a/., 1988, 1989). 

Replication dynamics becomes more complex if the restriction to frequency 
independent replication rate constants is dismissed (as we do here). Frequency 
dependent replication functions are readily interpreted in terms of chemical reac- 
tion kinetics: positive sign means catalysis, negative sign implies inhibition. In the 
simplest case the catalyst is a molecule of the same class as the template. Exam- 
ples are now well known in the biochemical literature: certain RNA molecules can 
act as specific catalysts for processing (Cech, 1986) or replication of other RNA 
molecules (Doudna and Szostak, 1989). Catalyzed replication is thus dealing with 
two forms of catalytic action (section 2): one RNA molecule is copied and acts 
as autocatalyst or template and a second one has the catalytic function (like a 
conventional protein enzyme in conventional biochemistry). Biochemical cataly- 
sis usually follows complex multistep mechanisms whose detailed dynamics is not 
accounted for by the single step kinetics of replicator equations. Replicator dy- 
namics, however, provides straightforward insight into phenomena observed under 
some limiting conditions. As an example we mention "faster-than-exponential" 
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growth predicted by the hypercycle equation (a special class of replicator equa- 
tions discussed in section 2; see also Eigen 1971, Eigen and Schuster, 1979) as 
recently found experimentally with replication of RNA viruses in host cells (Eigen 
et a/., 1991). 

Frequency dependent error-free replication has been studied extensively in the 
past (for a comprehensive survey see Hofbauer and Sigmund, 1988). All kinds of 
complex dynamical behaviour including chaotic attractors (Schnabl et a/., 1991) 
were observed. Introduction of mutation into these dynamical systems provides 
substantial difficulties for qualitative analysis (see section 2) and very few investi- 
gations with explicit and exact consideration of mutations were performed so far 
(for a general approach to mutation based on perturbation theory see Stadler and 
Schuster, 1992). There is one important exception: the selection mutation equa- 
tion derived from Fisher's selection equation by explicit consideration of mutations. 
It leads to much simpler dynamical phenomena as oscillations and deterministic 
chaos can be excluded, and analytical results are available on mutations in this 
particular case (Hadeler, 1981; Hofbauer, 1985). 

In this paper we shall derive exact results complemented by numerical studies 
on autocatalytic replication networks including explicit mutation terms. At first 
we introduce the replication-mutation equation in section 2 and consider the rest 
point in the central part of the physically meaningful domain of the phase space. In 
section 3 the two species model is presented in great detail. Section 4 and section 5 
are dealing with mutation in special cases of replicator equations, in the hypercycle 
model (Eigen and Schuster, 1979) and in the multidimensional generalized Schlogl 
model (Schlogl, 1972; Schuster, 1986). In the latter case we had to restrict the 
analysis to low dimensions (n = 3 and 4). A few results are available also for the 
limit of large dimensions (limn — > oo). 
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2. The replication-mutation equation 

2.1. The model and its assumptions 

The fundamental process basic to the simplest class of autocatalytic reaction net- 
works with frequency dependent replication rates is of the form 

(A) + I, + Ij- y Ik + I^ + IJ■, z,j,A; = l,...,n . (1) 

A is the substrate needed for replication. Its concentration is assumed to be 
buffered and does not enter the kinetic equations as a variable (and is hence 
written in parentheses). The replicator produced in the reaction is the target lyt, 
li is the template which is replicated, and Ij is the catalyst, all three belonging 
to the same class of species. The non-negative rate constants for these reactions 
are understood as elements of an n X n matrix A = [aij > 0). We assume that the 
species (the term species is used here for defined biochemical or chemical entities 
entering the kinetic equations as individual variables and not in the biological 
sense) 1^ replicates correctly with the frequency qa. A mutation from 1^ to lyt, 
[k ^ i) occurs with the frequency q^i- The product of a replication process has to 
be either a correct copy or a mutant and hence the conservation law 

n 

Y^qM = I (2) 

k=l 

holds. Thus, the mutation matrix Q = {qui) is a (column) stochastic matrix. 
Mutation frequencies depend on template I i and target lyt, but not on the catalyst 
Ij. The relative concentration of a species is denoted by [Iyt]/(7 = Xyt, where 
^ = YJk=i[^k\, therefore YJl^-^ Xk = 1. 

Application of mass action kinetics to the reaction-mutation network (l) under 
the constraint of constant total concentration C yields the replication-mutation 
equation, an ODE of the form 

n n 

Xk = "Y'^qkz a^j x^Xj - Xk ^ ; k = l,...,n. (3) 

1=1 j=i 
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An unspecific but time dependent dilution flux (Eigen and Schuster, 1979) 

n n 

Z=l J = l 

has been introduced for keeping the sum of relative concentrations normalized to 
unity at every instant. This constraint results in an internally controlled constant 
population size and is tantamount to a common assumption in population genetics. 
Since the rate constants aij are non-negative the flux is non-negative as 

well. This particular form of the flux has been found to be most convenient for 
mathematical analysis. In addition, the qualitative dynamics is essentially the 
same for a continuously stirred tank reactor (CSTR), where the dilution flux is 
constant (Schuster and Sigmund, 1985). 

By deflnition relative concentrations are non-negative quantities and hence 
equation (S) is physically meaningful on the simplex 



n 

Sn = ^{xi,. . . ,Xn) I a:^ > 0, y^x^ = l| 

i=l 



only. If aij > 0, for all j = 1, . . . n then Sn is a compact, forward invariant set 
for the ODE (3). 

If we choose the special case of error-free replication {Q = E, the identity 
matrix), the dynamical system (S) collapses to the second order replicator equa- 
tion^ 

n n n 

Xk — Xk ^ ^ ^ Oj^i Xl ^ ^ ^ ^ Ojlj X^X (4) 
i=l i=l j=l 

(for a review see Hofbauer and Sigmund, 1988) which we for convenience in vector 
notation, 

i = 7^(x) , (4') 

with ). The replicator Held 7^(x) = (T^i, . . . ,7^n) is deflned in 

terms of its components 

n n n 

7^yt(x) = Xki^'^aktXt -"^"^a^j x^Xj^ . (5) 

i=l i=l j=l 



P.F.StADLER et.al. : AUTOCATALYTIC NETWORKS 



Page 7 



In order to make the replication-mutation equation more handy we introduce 
a mean replication accuracy a mean error rate p, and the mean mutation rate 
e by 

1 " 

q = - S^ q^^ ^ p = 1 -q and 



n 
1=1 



n n 



^ ^ i=i ^ ^ i=i j=i 

According to the conservation law (2), p, q, and e can be converted into one 
another by the following expression: 

p = 1 — q = (n — l)e. 

The forthcoming analysis is facilitated by the definition of a modified mutation 

matrix W = (wij) whose off-diagonal elements are weighted against the mean 

mutation rate, and whose diagonal elements are set to zero: 

for z = J , 

q^J I e for I ^ J . 

The matrix W is used to define a mutation Held Ai{x.) = (A^i, ... , A^n) with 

n n 

A^yt(x) = "^"^{wkz a^j - w^k akj Xk) Xj . (6) 

1=1 j=i 

We note that the terms with i = k cancel and hence the assumption of vanishing 
diagonal elements, t/;ytyt = 0, has no influence on the mutation fleld. 
The model equation (S) takes now the form 

i = £(x) = 7^(x) + eA^(x) . (3') 

We are dealing with circulant reaction matrices A = {aij = ai^ij^i; z,jmod n} 
predominantly, and hence it is convenient to measure all reaction rates relative to 
the diagonal term which represents the autocatalytic rate constant an = a. This 
means that we split 

Oil = aij — a, with an = Q . 
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After dropping the hats, ciij =^ aij, we obtain from equation (3') 

i = 7^(x) + eA^(x) + eaB(x) , (7) 
where the background Held B{x.) is defined by 

n 

Byt(x) = "^{wkzX^ - W.kXk) ■ (8) 
i=l 

Remark. Equation (7) reduces to x = eaB(x) in the absence of selection, i.e., if 
A = 0. This mutation equation has a unique and globaUy stable interior equilib- 
rium if ae > and the mutation matrix Q is irreducible (Hofbauer and Sigmund 
1988, p. 251). 

Remark. If the mutation matrix Q is bistochastic then we can use the following 
shorter form for mutation and background fields: 

n n n 

Mk{x) = '^Wkz'^Xj {a^jX^ - GkjXk) and gfc(x) = y^^Wkzjxt - Xk) ■ 

i=l i=l i=l 

This includes the important special cases of both symmetric and circulant mutation 
frequencies. 



2.2. Existence of a central equilibrium 

Since Sn is compact and forward invariant for the replication-mutation equation, 
Brower's fixed point theorem ensures that there is always a rest point in the interior 
of the simplex. 

Many replicator equations (4) sustain an equilibrium in the interior of the 
simplex Sn- This fixed point can be readily shifted into the center c = ^(1,...,1) 
by a barycentric transformation (Hofbauer and Sigmund, 1988). It is useful to 
derive the conditions which leave the central fixed point unchanged. 
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Lemma 2.1. Let c = ^(1,...,1) be an equilibrium point of the replicator 
equation (4). Then c is also a fixed point of the replication-mutation equation (S) 
if the condition 

WAl = -(lAl) -Wl (9) 

n 

is fulfilled. 

Proof. Evaluation of equation (S) at the central fixed point, TZ{c) = 0, yields 

n n n 

i=l i=l j=l 

or in matrix form n[Al)k = (1^1) for all k. Substituting this into the replication- 
mutation equation (S) at c, 

_ n n n n 

Xk{c) = ■^(^'^'Wkza^j - y^^Wki y^gfcj), 
^ i=l j=l f=l j=l 

and rewriting in matrix form completes the proof. ^ 

Theorem 2.1. Suppose A and W are circulant. Then c is is a rest point of the 
selection mutation equation (S). 

Proof. Let C' = The product of two such matrices fulfills C'C^ = C""+' 

since 

^ \ki — / ^ Ok,j — lOj,i — m — Ok,i — m — l — 

i 

where indices are defined modulo n. If A and W are circulant they can be repre- 
sented by 

n—l n—1 

A=Y^ GjC^ and W = J2 "^jC^. 

j=o J=0 

Then WA = w.a^-.C"" and C'l = 1 implies 
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for arbitrary k. Since A is circulant we have (Al)i = ^(lAl), and the condition 
in lemma 2.1 is fulfiUed. ^ 

In order to create a fixed point at c it is by no means necessary to fulfil the 
condition of theorem 2. For instance, if the reaction matrix A is a multiple of the 
identity matrix E, then c is a rest point of equation (S) for arbitrary mutation 
matrices Q. 

In the forthcoming analysis we shall use different assumptions for the mutation 
matrix Q in order to be able to derive analytical expressions for the locations and 
the stability properties of rest points. These assumptions are essentially related 
to two basic models: uniform mutation rates and uniform error rates. 

2.3. The uniform mutation rate model 

The simplest and most restrictive, as well as most unrealistic, model for mutation 
assumes that all mutation rates are equal, Wij = 1 for i ^ j. We will refer to this 
mutation matrix as Q*^"^-*: 

qW = 

The corresponding W matrix is the of the form W*^"^-* =1 — E. Here I denotes the 
matrix with all entries equal to 1. 

Lemma 2.2. The central point c is an equilibrium of the replication-mutation 
equation (7) with mutation matrix Q if c is an equilibrium of the replicator 
equation (4). 

Proof. Substituting Q*^"^-* =1 — E into lemma 2.1 yields 

lAl -Al = (lAl)l -Al = (lAl)l - -(lAl)l, 

n 

i.e., we need {Al)k = ^{lAl) for all k. This condition is necessary and sufficient 
for c to be an equilibrium of the replicator equation (4). _ 



1 — [n — l)w for i = j 
w for i ^ j 



P.F.StADLER et.al. : AUTOCATALYTIC NETWORKS 



Page 11 



2.4. The uniform error rate model 

The model assume that a species I yt is represented by a sequence of length u. 
We further assume that mutants appear through copying one or more symbols 
incorrectly, and that the error rates are independent of the position of the digit 
in the sequence (Eigen, 1971; Eigen and Schuster, 1977; for a recent survey of the 
uniform error rate model see Eigen et a/., 1988 and 1989). The fidelity q of the 
replication process is a measure of the probability that a single symbol is copied 
correctly. The Hamming distance d[i^k) of two sequences 1^ and Ik counts the 
number of digits in which they differ (Hamming, 1950). Binary sequences of length 
u can be readily mapped on a hypercube such that the Hamming distance of the 
sequences counts the minimum number of edges that have to be passed in order to 
go from one sequence to the other. For short we call such a model the hypercube 
model and we obtain a specific form for the mutation matrix Q, which we shall 
denote by Q 

For convenience we introduce b = and the matrix B{u): bij = fe'^^*'^) for a 
i^-digit sequence. B[u) is a 2^ X 2^ matrix. The B-matrices can be obtained by 
recursion (Rumschitzky, 1987): 

The eigenvalues of B[l) are \± = 1 it 6 with the eigenvectors X-t = (1, ±1). The 
eigenvectors x_ and x_|_ are orthogonal. Let us assume that Xi is an eigenvector 
of B{u) with the eigenvalue Then y_|_ = (xi ® Xi) and y_ = (xi Xi) are 

eigenvectors of B[u + 1) and for the eigenvalues holds \[u + 1)-|- = (1 it b)\[u). 
Thus we can write down the eigenvalues of B[u) explicitly: 



x^{u) = {i+br-^{i-b)\ j = o,...,// 
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where Xj[u) is (p-fold degenerate. The above procedure produces a complete set 
of orthogonal eigenvectors of B[u). In particular, if we choose the +-sign in each 
recursion step we find that 1 is a simple eigenvector of B[u) for all dimensions. Of 
course B[u) and Q*^^-'(i^) have the same eigenvectors. The eigenvalues of (^^'^\u) 
are readily obtained: 

</>, = (2g-l)^ (10) 



2.5. The single point mutation model 

The hypercube model is simplified further by the assumption that error rates are 

so small that only single point mutations occur frequently enough to be considered. 

This is tantamount to assuming (p = 1 — q = (n — l)e): 

= [I -pf I -2p l-2w 

q{l —q) = (1 — p)p ^p^w 

(1 - g)2 = p2 ^ ^ 
Mutations are thus restricted to the interconversion of species with Hamming 

distance d[i^j) = 1. The recursion reads: 

We find, that the eigenvectors are the same for B[u) and B[u). For the eigenvalues 
we find 

A,(//) = l + (//-2j>, j = 0,...,// 

Again the eigenvalue \j is (p-fold degenerate. For the matrix Q*^^-*' = q" ■ B we 
have =q^{l-\-vh) = q^ -\-vq^~ {1 — q). In order to keep the sum of mutation 

probabilities normed to one, we define: 

q(3) ^ \ q(3)/ 

+ //g^-l(l - g) 

We can now easily calculate also the eigenvalues (^j of Q*^^-'(i^). 

^ , + (.-2,)(l-,) 

q + v(l-q) ^ ' 
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3. Models with two species 

3.1. The general case 

The (bistochastic) mutation matrix of a replication-mutation-system with two 
species is uniquely determined, 

Q= (^^~^ 1 -e) ' "^'^^ < e < 1. 
The most general reaction matrix is 

A = (^'^ with a,6, c,c/>0. 

The selection mutation equation is defined on the unit interval. We can eliminate 
the concentration of species I 2, 3^2, by setting We obtain 

i =[-a + b + c - d]x^ + [a-2b- c + 2d+ {-a + b- c + d)e]x^ + 

(12) 

+ [b - d + {-b + c - 2d)e]x + de = fix) 

The rest points of equation (l2) can be obtained by solving the cubic equation 
f{x) = 0. It turns out, however, that the resulting expressions are too complicated 
to be useful for further analytical work. 

The mutation free system has been analyzed completely by (Eigen & Schuster, 
1978). There are four generic phase portraits 

a > c A d > b Competition of 1 1 and I 2 (bistability). 

a > c A d < b Selection of 1 1, i.e., one equilibrium at xi = 1. 

a < c A d > b Selection of I 2, i.e., one equilibrium at 2:2 = 1. 

a < c A d < b Cooperation of 1 1 and I 2 . 
The phase diagram for the mutation free case is show in figure 1 (a) in terms of 
the parameters 



r 



I -\- b -\- c — d 



and 



s = a-\-b — c — d 
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Figure 1: Phase diagram of the two species model, (a) for £=0, (b) for small £>0. 



It is possible to generalize these results to small mutation rates using a perturbation 
approach (Stadler and Schuster, 1992), see figure 1 (b). 

If a = d = then x = and x = 1 are rest points of equation (3') for all e. 
It can be shown analogous to the proof of lemma 3.2 below that they are always 
unstable. In the following we will assume that (a, c/) ^ (0,0). We will use the 
following abbreviations for the coefficients of f{x): 

as = —a -\- b -\- c — d = r 

a2 = a — 2h — c -\- 2d -\- [—a + 6 — c + c/]e 

(13) 

ai = h — d + [—h + c — 2d]e 
ao = de 

Lemma 3.1. Let e > 0. Then there are only two possible generic phase portraits 
consisting of 

[A] a unique globally stable rest point, or 

[B] two locally stable rest points separated by an unstable one 
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The bifurcation between the two phase portraits is genericaUy a saddle-node bi- 
furcation. 

Proof. The unit interval is strictly forward invariant, hence there is at least one 
stable rest point. If it is degenerate, i.e. in the case it has a zero eigenvalue, 
then we have a pitchfork bifurcation. Since f{x) is cubic, there are at most three 
equilibria and at most one of them can be degenerate. If a rest point is degenerate, 
then there are at most two equilibria (this situation corresponds to a saddle-node 
bifurcation, or in a degenerate case to a transcritical bifurcation). If there are 
exactly two equilibria in (0, 1) then continuity ensures that the flow changes its 
direction at only one of them, hence the second rest point is degenerate. This 
situation is of course not generic since an arbitrarily small perturbation of any of 
the coeflicients leads to a phase portrait with either one or three rest points. If 
we have three rest points in (0, 1) then none of them can be degenerate and, by 
continuity, two of them are stable and the unstable equilibrium lies between them. 

I 

Lemma 3.2. If r is non-negative, then equation (l2) has phase portrait [A] for 
aU e > 0. 

Proof. If as > 0, then equation (l2) has either a single unstable rest point 
somewhere on the real line (/' > everywhere), or it has three rest points on 
the real line, two which are unstable (/ has then a local maximum and a local 
minimum. Lemma 3.1 guarantees that the flrst possibility can never occur for a 
selection mutation equation. For the second case we know (again from lemma 3.1) 
that there must be a single stable equilibrium in (0, 1). A degenerate rest point 
occurs if f{x) = for the local minimum or the local maximum of /. In both cases 
the second rest point must be unstable, contradicting lemma 3.1. For the special 
case as = we have at most two rest points on the real line. Suppose f{x) = 
has two real solutions, then one of them corresponds to a stable rest point, the 
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other one to an unstable one. Lemma 3.1 guarantees that we have phase portrait 
[A]. More degenerate cases cannot occur since we have always a stable rest point 
in (0,1). . 

Lemma 3.3. A necessary condition for the occurrence of phase portrait [B] is 

(i) — 3a3 > a2 > 0, and 

(14) 

(ii) ai + a2 > 0. 

Proof. For two stable equilibria existing in the unit interval it is necessary that 
f'{x) has its maximum in (0, 1) and that f'{xm) > 0, since there must be an 
unstable equilibrium in (0, 1). From 

f"{xm) = Qa3Xm + 2a2 = 

we find x^ = —02/303. Lemma 3.2 implies that 03 < is fulfilled. Then x^ £ 
(0, 1) implies condition (i). A short calculation shows that f'{xm) = ci2Xm + xi. 
Since < 02 < 02, it is necessary that at least 02 + oi > 0. ^ 

Corollary. There is an unique rest point for 

6 + c , , 

e>l —- 15 

~ a+d ^ ^ 

Proof. A sufficient condition for ruling out phase portrait [B] is ai + 02 = 
(a + (/) — (6 + c) — e(a + c/) < 0. Equation (l5) is an immediate consequence. ^ 

Hence there is a critical mutation rate e*(a,6, c,c/) < 1 such that there is a 
unique and globally stable equilibrium for all e > e*. Equation (l5) provides an 
upper bound on e*. This bound is reasonable, in the sense that it reproduces the 
condition 03 < which is necessary for the existence of phase portrait [B], but 
which is far from being optimal. 
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Remark. For e = 1/2 we find that x = 1/2 is a rest point for all coefficients 
a, 6, c, d. 

The bifurcation points of equation (l2) as a function of e can be found in 
principle by solving 

fix) = A fix) = 0. 

In case of the general model we obtain a quartic equation for e with no simple 
solutions. 



. B 



. 4 



. 2 



























1 / /{/ 
1 //// 

- 1 '/# 

^ 1 1 1 1 1 1 1 1 1 1 1 1 


1 1 1 1 1 1 1 1 1 1 1 


1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 



.2 ,4 .6 .SeI-D 



Figure 2: Numerical bifurcation diagram showing the solutions of f(x)=0 for a = 10,6=5, 

(1=2(10 — c), and 0<c< 10 as a function of £. Each choice of c is shown in its line style. The 
dashed parabola belongs to the symmetric case c=5, for which we observe a pitchfork 
bifurcation. The dot-dashed parabola £=1/4 — (1/2 — x)^ defines the boundary of the 
region within which we found more than one rest point in all our numerical examples. 

For small mutation rates we need a > c and d > b for the existence of three 
equilibria. Numerical solutions of f{x) = indicate that these conditions are 
necessary for arbitrary e. We have calculated the bifurcation diagrams of such 
systems as a function of e for a large number of choices for a, 6, c, d (see figure 2 
for an example). Our numerical experiments suggest the following 
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Conjecture. e*{a,b,c,d) < 1/4. 



3.2. Special cases 

The additional constraint a+6 = c-\-d ensures that c is a rest point, i.e, /(1/2) = 0, 
for aU e. A rather lengthy calculation shows that this rest point is stable if 



- 1^ 

£ > T 1 7 

4 V a + d 



Even for this simplified model we did not succeed in finding a better bound for e*. 
The special case of a symmetric reaction matrix 



A 



a b 
b a 



(16) 



contains both the hypercycle A = and the two-dimensional Schlogl model 
A = E as special cases. It can be analyzed in full detail. 

S H 




b/a 



Figure 3: The parameter space for the symmetric two species model (l6). Shown are different 
phase-portraits giving rise to qualitatively different dynamics for the parameters ^ and 
s (Note: at ■|- = 1 and £=0 every point in [0,1] is a rest point). S and H denote the 
Schlogl model and the hypercycle model, respectively. 
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Table 1. Eigenvalues of the Jacobian at the fixed points. 









Ao 


1 

2 


a(2e - 1) 


Ai 


^ -2ae 


(6 — a) + 4ae 



Lemma 3.4. The symmetric model, equation (l6), has an unique and globally 
stable equilibrium, i.e, phase portrait [A] if 

For e < e* we have phase portrait [B] . 
Proof. Equation (l2) reduces to 

X = 2(6 — a)x'^ — 3(6 — a)x'^ + (6 — 2ae — a)x + ae. (17) 

with the equilibria 

1 1 jl ae . 

xq = x± = - ± \ 18 

2' 2V4a-6 

The square root in equation (l8) is real and non-zero for e < (1 — 6/a)/4, and 
and x± lies in [0, 1], hence we have three fixed points. For e = e* the square root 
vanishes and x± = xq = 1/2, i.e., there is only one rest point which must be stable 
by lemma 1. For e > e* the equilibrium x = 1/2 is unique. ^ 

Lemma 3.4 shows that that we cannot find a better uniform bound for e* 
than the conjecture e* = 1/4. 

It is straightforward to compute explicitly the Jacobian and its eigenvalues 
at the positions of the fixed points (Table 1). The eigenvalue Aq belongs to the 
external eigenvector which does not influence the dynamical behaviour. 
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With e = e* and y = x — 1 /2 we find that equation (l7) reduces to 
The analytical solution of this ODE reads: 

1 

y{t) = sgn(y(0)) 



and (272) stable if a ^ b. If a = 6 we have e* = 0, i.e., we recover the mutation 
free replicator equation with all entries in the reaction matrix being equal. In this 
case all points are rest points. 

Figure 3 shows the regimes of phase portraits in the two-dimensional space 
of the parameters e and ^. For ^ — > we have the Schlogl model and for ^ — > 
00 we have the hypercycle as limiting case. Obviously this parametrization is 
sufficient to describe our model. The parameter ^ measures the cooperativity in 
the system: the Schlogl model refers to no cooperativity = O), and the hypercycle 
is "cooperativity pure" (lim - — > 00). 



4. The hypercycle model 

The reaction matrix A of the elementary hypercycle (Eigen and Schuster, 
1978) is given by , i.e., by 

a,j = S,^j-i; z,j = 1,2, . . . , modn , 

and represents a special case of a (non-symmetric) circulant matrix. The corre- 
sponding replicator equation has an asymptotically stable rest point in the interior 
of the simplex Sn for n = 2, 3, and 4. For n > 5 the central fixed point becomes 
unstable and there is an asymptotically stable limit cycle in the interior of Sn 
(Hofbauer et a/., 1991). 
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The replication field 7^(x) and the mutation field A^(x) in the hypercycle 
model simplify to 



7^yt(x) = Xk[xk-i - '^XiXi-ij and 

(19) 

Mk{^) = "^Wki {xiXi-i - XkXk-l) , 

t=l 

where all indices are understood modn. The background field -B(x) remains un- 
changed. By theorem 2.1, c is an equilibrium of the replication-mutation equation 
of elementary hypercycles for arbitrary n, a, e and W. 

As shown in (Stadler and Schuster, 1992) there are no other rest points for 
small e > 0. Numerical studies indicate that the hypercycle model with mutation 
has in fact a unique globally stable equilibrium for all e > 0. In the following we 
will restrict our attention the the rest point c. 

The Jacobian at the central equilibrium, D = [dij = dCi/dxj[c) ^ ^ is readily 
computed. Its elements dki are of the form 

dkfic) = - (Sk-1,£ - -) + - (wk£ + Wk,£+1 - (n - l){Sk,£ + Sk-1,£)) + 

n \ n / n \ / ^20) 



+ ae(^Wk£ - [n - l)Sk,£^ . 



The properties of the central fixed point of elementary hypercycles with mutation 
are summarized in several lemmas and theorems. 

Lemma 4.1. D is circulant if and only if W is circulant. 

Proof. Define Yki = {a + -)wu + -Wki+i and 

n n 

1 2 1 — n 

Zkl = ) H £{h,l + - Oi£{n - 

n n n 

We have D = Z -\-'eY . Z is circulant since it is a linear combination of the matrices 
E = (Sij) and = {5i-ij) which are both circulant. Clearly, D is circulant if 



P.F.StADLER et.al. : AUTOCATALYTIC NETWORKS 



Page 22 



and only if Y is circulant. It is obvious that Y is circulant if W is circulant. It 
remains to show the converse, i.e., suppose Y is circulant. Then 

Yu = (a + -)wn + -wii+i = (a + -)wi+ji+j + -wi+ji+i+j = Yi+ji+j 
n n n n 

must hold. Setting 1 = 1 and using wn = we find t/;i2 = 'Wi-\.j^2-\-j for all j . Now 
we proceed by induction; suppose wi^i = wis^j^is^j. W is circulant if and only if 

r , r 1 1 
(a + -)wi^i - (a + -)wij^j^ij^j = -wij+i wi+jj+i+j = 0. 

Tlj ft ft ft 

This simplifies to wi^is^i = i.e., W must be circulant. ^ 

Lemma 4.2. If W is circulant then the real parts 7^ of the eigenvalues (f)^ of 

D for m = 1 , . . . , n — 1 are given by 

1 27rm _n — 1 27rm _ 

7m = — cos e (1 + cos ) — aein — 1) 

n n n n 

e , 27rm(j - 1) _ ^ 27rm(j - 1) (^l) 

H > , + ,+1) cos h ea > w^jcos 

n ^-^ n ^-^ n 

1=1 j=2 

The remaining eigenvalue (j)o = —1/n belongs to the eigenvector 1 and does not 
influence the dynamics. 

Proof. Let cj = di^k,j+k- If D is an arbitrary circulant matrix, the eigenvalues 
of D are given by (f)m = Sj=i '^j^ri^'' where A„ = exp(^z). Hence 



Sftcf>^ = J2 Cj^Xn^'~'^ = J] C, COS 



27rm(j — 1) 



n 

1=1 1=1 



A rather long but straight forward calculation then yields equation (2l). In order 
to check that 1 is an eigenvector we calculate 



-{h-i,l ) + -{wki + Wk,i+i - [n - l){Sk,i + Sk-i,i)) + 

n n n 

1=1 

+ ae{wki - [n - l)5k,i) = 

1 

n 
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This eigenvector 1 points outside the simplex. Since Sn is invariant for equa- 
tion (4), the dynamical behaviour of our system is determined by the remaining 
n — 1 eigenvectors , . . . , (j)n-i ■ | 

Theorem 4.1. If the mutation matrix W is circulant, the central fixed point of 
the 2-, 3- and 4-species hypercycle is stable for a,e > 0. In the case n = 4, a = 
we need the additional condition W12 7^ 3, i.e. more then one mutation probability 
Qij} ^ 7^ i must be nonzero. 
Proof. 

(1) n = 2 The dynamical behaviour is determined by the eigenvalue 

71 = -- - -W12 - 2ae 

Obviously 71 < holds under the conditions of the theorem, therefore 
c is stable. 

(2) n = 3 We get one double eigenvalue 

1 e,3 , 
71 = 72 = -g + gig"'!^ - 3) 

Since W12 ^ 2 the bracket is nonpositive under our conditions and we 
find 71 < 0. 

(3) n = 4 In this case we get one single and one double eigenvalue 

71 = 73 = -2(^^12 - 3) - ae(3 + W13) 
72 = - 2ae(3 - wis) 

Using Wij < 3 we have 72 < and with W12 < 3 71 also remains 
negative for e > 0. For the non-mutating four membered hypercycle 
(e = 0) we know that the inner fixed point is stable (Eigen and Schuster, 
1978; Schuster et a/., 1978).. 
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For larger hypercycles the expressions get much more comphcated due to the 
increasing number of mutation terms. We therefore restrict our investigation to 
the special mutation matrix Q*^"^-*. 

Theorem 4.2. Let Wij = 1 for i ^ j and n > 5. Then bifurcations at the central 
fixed point c occur at e = £^(0; n-) with 



1 cos ^ n 



n 1 + cos + na 



1 < m < 



4 



{22] 



The central equilibrium c is stable if e > ei(a; n). If n = (mod 4) then there is 
a bifurcation for e = 0. 



e = 0.03 



£ = 0.05 




Si. 0.20 




0.20 

x(1) 



Figure 4: Bifurcation of the hypercycle with n=5. The limit cycles on the Ihs corresponds to 
error rates £ below the critical value. The trajectory on the rhs represents a typical 
solution curve above the critical value. It spirals inwards and ultimately reaches the 
stable rest point in the center of 5*5. 



Proof. Because of the simple form of Q expression (22) simplifies considerably. 
Substituting the identities 

27rm(j — 1) ^^-4 27rm(j — 1) 27rm 

cos = — cos 



n = ^ 

J = 2 

into equation (22) yields 



n — ' n n 

J=2 J=l 



1 27rm _ 27rm 

7m = — cos e[l + cos h Q,n\ 

n n n 
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The expression in brackets is always non-negative. Thus 7^ is always negative if 
cos < 0. Therefore we find simple zeros of 7^ for l<m<jate = £^(0; n). 
For 1 < m < n — 1 we have cos ^ > cos therefore ei(a; n) > £^(0; holds 

for 1 < m < J. Since c is stable if all eigenvalues 7^ are negative, we find that 
e > ei(a; n) is sufficient condition for stability. 
If n = mod 4 then e„/4(a; n) = 0. ^ 

We remark that and (j)n-m form a pair of conjugate complex eigenvalues 
if and only if m ^ -j: (f)jn = 7m + ^ ■ ^^^^^ sin Therefore the bifurcations at 
£^(0; with 1 < m < J are Hopf bifurcations. Numerical work for n = 5 shows 
that the bifurcation at e = ei(a; 5) ~ 0.047213 is supercritical figure 4. 

We mention further that a rigorous proof for the existence of asypmtotically 
stable limit cycles in the elementary hypercycle was already hard to derive for the 
error-free system (Hofbauer et a/., 1991). 

5. The generalized Schlogl model 

In Schlogl's model the replication matrix is defined to have only diagonal 
terms aij = • Sij. We consider here the simpler case of equal rate constants 
{ttij = Sij). From Lemma 2.1 follows that c is a fixed point of equation (7) for 
arbitrary a,e and W. Because of the simple form of A, it is more convenient to 
use equation (7) than equation (S). We have 

n n n 

i/t = ^ qkzX^t - ^k'^xj - qki{xr - xj,) (23) 

i=l i=l i=l 

The case of two species has been treated already in section 3.2. 
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5.1. Models with three species 

The calculations for the general three species case are rather sophisticated and 
thus we can treat here only the most simple model 

/100\ /l-2e e e\ 

A= 010 , Q = { e l-2e e 

\001/ \ e e l-2ej 

In this case the coordinates of all fixed points can be calculated for arbitrary e. 
From symmetry reasons, we expect three types of fixed points: 

(i) a unique inner fixed point |l 

(ii) three edge-type fixed points y_) 

(iii) three corner-type fixed-points y_|_) 

The analytical results, which were obtained by using of the software package 
REDUCE, are compiled in table 2. We find three types of phase portraits in this 
system. For < e < ^ we have a smooth deformation of phase portrait of the 
Schlogl model, for e > 3 — 2\/2 there is a unique interior sink. Between ^ and 
3 - 2V2 ^ 0.17157 we find a rather strange phase portrait with four coexisting 
sinks, and three very narrow channels connecting small zones at the boundary 
with the central sink (figure 5). 

We find that the sinks at the corners and the saddles on the edges move in- 
wards with incerasing error rate. The saddles collide with the interior source at 
e = ^. As e increases further, the interior fixed point is left as sink, which is sepa- 
rated by the three edge-type saddles from the other three corner-type sinks. When 
e reaches 3 — 2\/2, the saddles collide with the corner-type sinks and annihilate. 
If e increases further, we are left with a unique stable fixed point. We illustrate 
this behaviour in a constraint-response plot in (figure 6), 
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Figure 5: Fixed points and selected trajectories for the tree-species Schlogl model. All trajec- 
tories start at the boundary of the simplex. 




Figure 6: Bifurcation diagram for the rest points of the model (J is the interior equilibrium, C 
denotes rest points evolving out of the corners of the simplex, and E labels fixed points 
starting at the edge of the triangle). Solid line denote sinks, dashed lines mark saddles, 
and dotted lines stand for sources. 
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Table 2. Fixed points and signs of the eigenvalues for the three-species Schlogl 
model 



Type 


Interior 


Edge 


Corner 


Coordinates 


|1 


{x+,x+,y-) 
x+ = 1(1 + e + p) 

y- = 1(1 - e - p) 


{x-,x-,y+) 
X- = 1(1 +e -p) 

y+ = 1(1 -e + p) 


Def. 


[0,1] 


[0, 3 - 2^2] 


[0, 3 - 2V2] 


Ao 


1 

3 


(l-3£)p+/ji 
4 


(3£-l)p+/ji 
4 


Ai 


l-2e 


(l-3?)p+/j2 
4 


(3£-l)p+/j2 

4 


A2 


= Ai 


3(l-3£)p+/j3 

4 


3(3£-l)p+/j3 

4 



p = \/e^ — 6e + 1 

/zi = -3e2 + lOe - 3, = -3e2 + 18e - 3, /zg = -Qe^ + 6e - 1 
5.2. Models with four species 

As in the previous section we choose the identity matrix as reaction matrix. We 
will investigate three different mutation models. The first one describes a one digit 
string with four possible characters, the second one corresponds to a binary two 
digit chain, and the third one corresponds to a binary chain, where mutations with 
Hamming distance d[i^j) > 1 are forbidden. 

5.2.1. Uniform mutation rate model 

In the first case all mutation rates are equal and the mutation matrix reads 
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Figure 7: Fixed points and some trajectories for the Schlogl model with four species and the 
mutation matrix VF=Q'^^-' (uniform mutation rate model). 




Figure 8: Bifurcation diagram for the Schlogl model with four species and the mutation matrix 
Q=Q'^^-' (uniform mutation rate model). Notation as in figure 6, and S denotes rest 
points developing from equilibria on the surface of the simplex. 
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Table 3. Fixed points and eigenvalues for the Schlogl model with four species 
and the mutation matrix given by Q = Q (uniform mutation rate 
model). 

Type Interior Edge Corner Surface 



Coordinates jl 



X 



y 



4 



, ^ l + 2£-2g2 

6 

, ^ l-2£+2g2 

2 



l + 2£+2g2 

6 

l-2£-2g2 

2 



Def. 


[0,1] 


[0,1] 


[0,1 


2 J 


[0, 1 - % 


Ao 
Ai 

A2 


1 


2e-\ 
1-8? 


2(4?- 


■l)?2 + fei 


2(l-4?)g2 + /ii 


4 

1-8? 


2(4?- 


3 

-i)g2+'*2 


3 

2(l-4?)g2 + /i2 


4 

= Ai 


2 

4£-l^ 

2 ?1 


4(4?- 


3 

-i)g2+'*3 
3 


3 

4(l-4?)g2 + '*3 
3 


A3 


= Ai 


= -A2 




= A2 


= A2 



qi = VI - 8e, ^2 = - 2e + I 

/zi = -8e^ + lOe - 2, /z2 = -8e^ + 16e - 2, /z3 = -16e^ + 8e - 1 



In the mutation field this matrix is tantamount to the mutation matrix W*^"^-*. 
Because of the high symmetry of our problem, the calculation of the fixed points 
simplifies considerably. There are four classes of fixed points (table 3). We have 

1 central fixed point (1;, 1;, 1;, 1;) 

4 corner type fixed points (x, x, x, y) 

6 edge type fixed points (x, x, y, y) 

4 surface type fixed points ( 
with all permutations of x^y being allowed. 

We find that bifurcations occur for e = i = 0.125 and for e = 1 — 
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0.133975. For e = the interior fixed point is a source, the corner-type equihbria 
are sinks, aU other fixed points are saddles. When e increases, all fixed points move 



points collide with the central equilibrium. The edge type points disappear, and 
one eigenvalue of the surface-type points chances sign, but the equilibria remain 
saddles when the come out on the other side of the interior fixed point. As for 
the three species model, the interior equilibrium becomes stable. We have five 
coexisting sinks, four of which are the corner type. Their basins are separated 
from each other by a narrow zone with tetrahedral symmetry that form the basin 
of central equilibrium. For e = 1 — the corner-type and the surface-type fixed 
points crash and vanish. Figure 7 gives some phase portrait for various values of 
e. In figure 8 we give a constraint response plot for this model. 

5.2.2. The hypercube model 

Individual species in this model are indentified with binary strings, for example we 
have 1 1= (00), I 2= (01), I 3= (10) and 14= (11) for n = 4, and the uniform error 
rate assumption is applied (since we 2^ binary strings of chain length u hypercube 
mutation models exist only for dimensions that are powers of two: n = 2, 4, 8, . . .). 

The hypercube model looks more involved at first glance. Nevertheless, we 
obtained fairly simple expressions for the coordinates and the eigenvalues of all 
fixed points. The we write down our mutation matrix, explicitly written as 



is circulant and we must have cyclic symmetry in the coordinates. In this model, 
however, the edges of the simplex 5*4 are no longer equivalent. Instead, there are 



the remaining two edges connect corners with d[i^j) = 2. We refer to them as 



into the interior of Sn- At e 



g the six edge-type fixed and the surface type fixed 




two classes: four edges connect 
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(01) 




(00) 

A 




(11) 



(10) 




Figure 9: The hypercube mutation rate model for four species. Full arrows (< >) refer to 

single point mutations, and dotted arrows (< • • • >) indicate simultaneous changes of 
two digits at a time. 



direct and indirect edges, respectively. There are also two different edge-type fixed 
points, and thus we are dealing with five classes of fixed points (table 4). 

1 central fixed point (|, |, |, |) 
4 corner type fixed points (x, y, x, z) 
4 direct edge type fixed points (x, x, y, y) 

2 indirect edge type fixed points (x, y, x, y) 
4 surface type fixed points 

For the mutation matrix Q ^'^^ there exist four different bifurcation points: 

51 = ^^^ 0.866, ^2 = 1(1 + ^) ~ 0.8536 

93 = f, ^4 = 1(1-^) - 0.1464. 

For q = 1 we have the unperturbed Schlogl model. As q decreases, all fixed points 
move inwards. The central equilibrium is a source, the corner type fixed points 
are sinks, and we have 10 saddle points in 5"^. At q = qi two pairs of surface 
type saddles collide with indirect edge type fixed points, one eigenvalue of which 
changes sign, and disappear. The indirect edge type equilibria remain saddles 
until q = q2 where they vanish crashing into the central fixed point. At this 
critical value of q the central equilibrium becomes a saddle, and remains so till 
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X4 X4 X4 




q=0.aODOO q=0.50QDQ q-0. 10000 

Figure 10: Fixed points and trajectories for the four species Schlogl model with mutation matrix 
Q=Q (hypercube model). 
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Figure 11: Bifurcation diagram for the four species Schlogl model with mutation matrix Q=Q 
(hypercube model). 
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Table 4. Fixed points and eigenvalues for the four-species Schlogl model with the 
mutation matrix given by Q^^^^^^ (hypercube model). 



Type Interior 


Corner 


direct 


indirect 


Surface 






Edge 


Edge 





Coord. il {x,y,x,z) {x,x,y,y) {x,y,x,y) {x,y,x,z) 

1 l + Vl l + '/2 Q 

A ±yA ^ A ^ A 2g-|-l 

_ _ _ l-g2 _ l + g3 

y — 2 y ~ 4 y ~ 4 y ~ 2(2g+i) 

r — r — l~g3 

^ ~ 2 ^ ~ 2(2g+l) 

Def. [0,1] [f,l] [f,l] [0,i-^]U [f,l] 



Ao 


1 

4 




/ 
2 




f 
2 


/ 

2g+l 




Ai 


4g-3 
4 


-f{4q - 3) 


4g-3 
2 




-8q+l 
2 


2q+l 




A2 


= Ai 


= Ai 


(4?-3)/ 
2 






2g^-g^-3g+2 + g4 
2g+l 


> 


A3 




1 -2g2 




-iif- 




2g^-g^-3g+2-g4 


< 


4 


2 


2g+l 



/ = 2g - 1, qi= v^^^, (Z2 = -8g+ 1, qs = 

qi = v^4g6 + I2g5 _ 19^4 _ 5^3 _^ 15^2 _ 6g + 1 



q = q3. Here the corner type and the direct edge type fixed points disappear by 
colliding with the central equilibrium which becomes stable. Between qs and q^ 
the central equilibrium is stable and unique. At q = q^ it loses stability again and 
becomes a saddle; two indirect edge type stable fixed points split off the interior 
equilibrium and stay stable until q = 0. We illustrate this behaviour in figure 9, a 
constraint response plot for the binary chain model is given in figure 10. 
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5.2.3. The single point mutation model 

In the third and last model mutations are restricted to single point mutations. It 
is a good approximaion for small error rates. The mutation matrix then reads 

/I — 2w w w \ 



w 1 — 2w w 
w 1 — 2w w 

w w 1 — 2w / 



where we use the parameter w = instead of e for convenience of calculation. 
The behaviour is very similar to the previous system, with the exception, that the 
central equilibrium does not become unstable for very low replication accuracy. 
As for the previous model, matrix Q ^'^^ we expect and find five classes of fixed 
points: 

1 central fixed point (|, |, |, |) 
4 corner type fixed points (x, y, x, z) 
4 direct edge type fixed points (x, x, y, y) 

2 indirect edge type fixed points (x, y, x, y) 
4 surface type fixed points ( 

The results for the single point mutation model are summarized in table 5. 

For w = we start with the same situation as above. As the mutation 
increases, all fixed points move inwards. The first bifurcation occurs at 

ao = 0.114078 

12^- VSTTl - 17^3 

where the surface type saddles collide with the indirect edge type saddles and 
disappear. We remark, that ao is the unique real zero of 

32x^ - 32x^ + 12x - 1 = 

The indirect edge type saddles change the sign of one eigenvalue and remain sad- 
dles; the collide with the central source at w = ^, disappearing and leaving the 
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X4 X4 X4 




p=0. 12500 p=0. 15000 p-0. 25000 

Figure 12: Equilibria and trajectories for the four species Schlogl model with mutation matrix 
Q=Q (single point mutation model). 




P-3/2E 



Figure 13: Bifurcation diagram for the four species Schlogl model with mutation matrix Q=Q 
(single point mutation model). 
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Table 5. Fixed points and eigenvalues for the four-species Schlogl model with 
mutation matrix q=q(^) (single point mutation model). 



Type 


Interior 


Corner 


direct 
Edge 


indirect 
Edge 




Surface 


Coord. 


■ ii 


(x, y, X, z'J 


{x,x,y,y) 


{x,y,x,y) 




(x, y, X, z'J 






_ 4«,--i+y7 

2(aw — S) 


X — ^'^'^^ 
4 


^ _ l-|-?2 

4 










_ v-\/f+\/9+h^/f 
y ~ 2(8w-3) 








y ~ 2(8w-3) 






^ 2{8w-3) 
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qi = VI - 4:W, q2 = VI - = -32w^ + 33w^ - lOw + 1 

/ = IQw* - Q4:w^ + 48w^ - 12w + 1, h= 16w^ - 12w + 4 
g = QAw'^ - 176w^ + 144w^ - 48w + 5, v = -4w^ + 8w - 2 
ao ^ 0.1140777468 



interior equilibrium as saddle point. For w = j the direct edge type saddles and 
the corner type sinks collide with the interior equilibrium and annihilate. The 
central fixed point is left unique and stable. We find this model with mutation 
matrix Q to have only one qualitative difference with the model with the Q 
matrix: The central equilibrium remains stable till the end of physical domain, 
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and the indirect edge type point do not reappear for large values of w (or e). We 
show some phase portrait in figure 11 and give a constraint response plot for this 
model in figure 12. 



0.2 0.4 0.6 0.8 1.0 

1 I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I C| 




Figure 14: Comparison of the stability of the fixed points for difli"erent mutation matrices. All 
mutation parameters have been converted to the mean replication accuracy f (Q 
uniform error rates; Q hypercube model; Q single point mutations). 



5.2.4. Comparison of the three models 

In order to compare the three mutation models the coordinates and eigenvalues 
of all fixed points were expressed as functions of the same parameter, the mean 
replication accuracy q = Qkk'- 

(1) uniform mutation rate model: q = 1 — (n — l)e 

(2) hypercube model: q = q'^ = q^'^^ 

(3) single point mutation model: q = q'^ = l — u-w = l — (Id n) ■ w 

The results are shown in figure 14. In all three models we observe a bifurcation 
at which a stable rest point jumps from a corner (representing the species with 
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the highest fitness value) towards the center. This bifurcation is tantamount to 
an error threshold of replication which will be discussed in the final section 6. We 
observe two interesting differences: 

(1) In the uniform mutation rate model this bifurcation is subcritical, wherease 
it is supercritical in the other two cases. Thus the uniform mutation rate 
model can lead to hysteresis phenomena. 

(2) In the hypercube model the rest point in the center becomes again unstable 
at very small values of q. Then the stable fixed point jumps onto the indirect 
edge. It corresponds to inaccurate complemetary replication (Swetina and 
Schuster, 1982; Schuster and Swetina, 1988; Stadler, 1991). In the limit 
g — > we have precise complementary replication: A + I®— ^Iq+I® and 

Apart from the stable rest point on the indirect edge (hypercube model) all rest 
points on edges and surfaces are unstable. 



D = i 

n 



5.3. The central equilibrium in high-dimensional systems 

The Jacobian for model (23) at the inner fixed point c reads 

2 + na)Q - (1 + na)E - -I . (24) 

n J 

Obviously, D is circulant if and only if Q is circulant. 

Theorem 5.1. The central equilibrium c of Schlogl's model, equation (23), with 
mutation matrix Q is a sink for e > ec* and a source for e < ec, where 

1 1 



n 2 + na 

Proof. W is circulant and qik = £ for i ^ k and qkk = 1 — {n — l)e for diagonal 
elements of Q. Hence equation (24) becomes 



dkl 
dkk 



[(2 + 



na e 



1](2 



for k ^ I and 
2- 



na e 



n 
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D is circulant due to Lemma 5.1 and therefore has the eigenvalues 

1^0 = dii + (n - l)di2 = , 

n 

jji = dii — di2 = —[1 — n{2 + na)e]. 

n 

It is easy to verify that D ■ 1 = —^1, i.e., /io belongs to the unphysical external 
direction, and the dynamical behaviour is determined by the (n — l)-fold degenerate 
eigenvalue /ii alone. ^ 

Lemma 5.2. Let Q be symmetric. Then each eigenvector of Q is also an 
eigenvector of D. The eigenvalues of D are 

IJk = - [(2 + na)(f)k - (1 + na) - 2Sk,o] 
n 

where (f)k are the eigenvalues of Q. 

Proof. Since Q is stochastic and symmetric 1 is an eigenvector, and all other 
eigenvectors Zyt can be chosen to be orthogonal to 1. Let I denote the matrix 
with all entries equal to 1. Multiplying the Jacobian, equation (24), with the 
eigenvectors Zk of Q yields 

1 2 
D -1 = - \(2 + na)P • 1 - (1 + na) • 1 1 • ll = 

n n 

1 2 
= — [(2 + na)cf)Q — (1 + na) • nl • 1 = /io ■ 1 

n n 

1 2 
D ■ Zk = - \{2 + na)P • z^, - (1 + na) ■ z^ 1 • z^,] = 

n n 

= - [(2 + na)(/)yt - (1 + na) - Ol • Zyt = /iyt • Zyt 

n 

I 

We can now determine the stability of the central rest point c of Schlogl's 
model with mutation terms defined by the uniform error model for arbitrary chain 
length V . 



Theorem 5.2. For the Schlogl model with mutation matrix Q hold 



s: 
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(i) The central fixed point c is a source if 



^ 2^ V 2 1 + 2^-ia^ 



(ii) The central fixed point c is a sink if 



1 . Jl 1 + 2^a . 1 1 1 + 2^a , 

2 - V 2 1 + 2^-1 J < ^ < 2^' + 2 1 + 2^-1 J- 

For the Schlogl model with mutation matrix Q holds: 

(i) The central fixed point c is a source if 

1 

q > 1 = Qu- 

(ii) The central fixed point c is a sink if 

2(2 + 2^a)-// 

< q < ——■ and 

2(2 + 2^a) - // + 1 

iy-4 

« > ^ I -I 
2^+1 

The second condition is always fulfilled if u < 4. 

Proof. The eigenvalues for the mutation models Q and Q*^^-*, respectively, are 
easily calculated from equation (24). We find 

n 

^M) = ^ • [(2 + 2^«)(29 - 1)^' - (1 + 2^«)] (25) 

- ( \ ^ + (;/ - 2j)(l - q) . 

/^.M = ^-[(2 + 2 a) -(1+2 a)] 

The eigenvalues iij[u)^ p,j become zero for 



=2('^V2-TT^' 

2j(2 + 2^a) - u 
~ 2j(2 + 2^a) -{y -I) 
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^[N /|\ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ? C 



10 12 

log2(n) 



16 



20 



Figure 15: Critical mean replication accuracy per replication round for different models and 
mutation matrices (* for mutation matrix Q O for Q and A for Q '^^-') depending 
on the number of species. Solid lines denote Schlogel's Model and the dashed line stands 
for the Hypercycle. 



where q- is a zero of /Uj(i^) if and only if j is even and qj is a zero in the physicaUy 
meaningful domain, i.e. qj G [0, 1], if and only if j > 2(2+2'^ a) • I 

The parameter q ist related to the mean replication accuracy q per replication 
round. For the full model Q*^^-*, the simplified model Q*^^-*, and the mutation 
matrix Q we find 



(2)1 



(3)1 



q + u{l- q) 



and 



(1)1 



in — l)e 



respectively. We are interested in the first bifurcation at the interior fixed point 
when we switch on the mutation, qc. Since we have used different mutation param- 
eters for convenience of calculation we transform them to q in order to compare 
them. For simplicity we set a = and n = 2" since the matrices Q and Q 
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can be defined only for powers of 2. We find 



2^+1 ' 

I' / -'- \i f 



For the hypercycle model with Q we find 

<;4<4 ]-l (1 z-'l +008 1(1-22-) 

We find that the parameter range in which the interior fixed point of the 
mutating system exhibits the same dynamical behaviour as in the pure replicator 
equation depends strongly on both the reaction matrix A and the mutation matrix 
Q. Note that the three different mutation matrices reduce to the unique matrix 
Q*^^-'(l) for n = 2, i.e, for u = 1. The values of Qc are compares in figure 14. 



6. Conclusions 

The results obtained for "catalyzed autocatalysis" (second order autocataly- 
sis: equ.l) with different model mutation models share the existence of an error 
threshold phenomenon with the simple replication-mutation system with frequen- 
cy independent replication rates (first order autocatalysis: Eigen, 1971; Eigen and 
Schuster, 1977; Eigen et a/., 1988 and 1989). In this simple case the system con- 
verges towards a unique and asymptotically stable stationary state provided there 
is no degeneracy of fitness values (i.e. the fitness of the fittest species called mas- 
ter species is larger than that of any other species) and the mutation matrix Q is 
positive definite. The mutant distribution is obtained as the largest eigenvector (i 
of the matrix Q ■ [a ■ E) where a = (ai, a2, . . . , a^), the vector of replication rate 
constants. From Frobenius theorem follows that all components of Ci are strictly 



P.F.StADLER et.al. : AUTOCATALYTIC NETWORKS 



Page 44 



positive. The eigenvector (i when computed as a function of the error rate p shows 
a threshold phenomenon at some critical p- value which is reminiscent of coopera- 
tive transitions. At low values of p the mutant distribution is centered around a 
master species (or master sequence if polynucleotides are considered). In terms of 
phase portraits this means that the unique stable fixed point P is located near a 
corner of the simplex Sn- As a function of p the equilibrium point moves slowly 
towards the center c. In the transition region around the critical error rate which 
is tantamount to the error threshold P jumps suddenly towards c and reaches the 
center asymptotically when replication becomes completely random (p=| in the 
hypecube model): lim-p^i/2 -P(p) = c. The longer the sequence {u), or the larger 
the number of species (n = 2^) is, the narrower becomes the transition zone, the 
sharper is the threshold. Depending on the vector of replication rate constants a 
the width of the transition region may remain finite in the limit of large numbers 
of species (limn — > oo) which is characteristic for a cooperative transition, or it 
may shrink to a transition point. Then we are dealing with a (higher order) phase 
transition (Leuthausser, 1987). 

The multi-dimensional Schlogl model (Schuster, 1986) is the higher order ana- 
logue to the simple selection case (with frequency independent replication rates). 
Instead of a unique stable stationary state we find n equlibria, one at each corner 
of the simplex Sn- Rate constants determine the sizes of the basins of attraction. 
Mutation causes the stable equilibria to move into the interior of Sn- At some 
critical error rate we observe a bifurcation (or several bifurcations) which make 
the central rest point to the unique stable steady state of the system. This is 
straightforward to visualize since the central rest point is the unique fixed point of 
the mutation matrix. Here, the error threshold resembles a first order phase tran- 
sition. This result was obtained for all mutation models, although some differences 
are observed in detail: in the uniform mutation rate model the bifurcation is sub- 
critical and the system shows hysteresis accordingly, whereas we find supercritical 
single jumps towards the center in the other two cases. 
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In general we observe that the onset of mutation simphfies rephcator dynam- 
ics: muhiple stable equilibria in the Schlogl model change into a single stable 
restpoint at sufficiently high error rates, oscillations in the mutation-free hypercy- 
cle model are replaced by the central equilibrium, and the chaotic attractor of a 
replicator equation with four species is turned stepwise into simpler dynamics if 
mutation terms are introduced. 
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